#!/bin/bash

set -o pipefail
set -o errexit
export LC_ALL=C

export BASE_DATA_PATH=/home/mixagol/data/
export BLAST_DB_PATH="${BASE_DATA_PATH}/1_databases/blast_db/GENOMES_DB"


export GENOMES_NUM=`cat "${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt" | wc -l`
g++ -O3 6_annotation/intersect_genes.cpp -o 6_annotation/intersect_genes
g++ -O3 6_annotation/intersect_genes_avg.cpp -o 6_annotation/intersect_genes_avg

export ANNOTATE_DIR="${BASE_DATA_PATH}/6_annotation/annot_$(date +%s)"
mkdir -p $ANNOTATE_DIR

cat gene.faa > $ANNOTATE_DIR/gene.faa

cat $ANNOTATE_DIR/genes.faa | 2_raw_matrix/run_blast.sh $BLAST_DB_PATH > $ANNOTATE_DIR/blast_results.txt

#for E_VALUE in 0.05 0.01 0.001 0.0001 0.00001 0.000001; do
for E_VALUE in 0.05; do

    export E_DIR=$ANNOTATE_DIR/e_$E_VALUE
    mkdir -p $E_DIR

    cat $ANNOTATE_DIR/blast_results.txt | awk -F'\t' '{if ($11<='$E_VALUE') print $1"\t"$2}' \
        | 3_normed_matrix/apply_maps.py \
            --genomes-map ${BASE_DATA_PATH}/3_normed_matrix/map_int_genom.txt \
        | 3_normed_matrix/create_matrix.py \
        > $E_DIR/vectors.txt
    
    6_annotation/intersect_genes \
        -m ${BASE_DATA_PATH}/3_normed_matrix/matrix_${E_VALUE}_.txt \
        -v $E_DIR/vectors.txt \
        -g $GENOMES_NUM \
        > $E_DIR/joined_raw.txt
    
    join -t$'\t' -1 1 -2 2 \
        <(sort -t$'\t' -k1,1 -S100M ${BASE_DATA_PATH}/3_normed_matrix/map_int_gene.txt) \
        <(sort -t$'\t' -k2,2 -S100M ${E_DIR}/joined_raw.txt) \
        | cut -f2,3,4 \
        > ${E_DIR}/joined.txt
    
    join -t$'\t' \
        <(cat ${E_DIR}/joined.txt | sort -t$'\t' -k1,1 -S100M) \
        ${BASE_DATA_PATH}/5_pathways/gene_kegg_map.txt \
        > ${E_DIR}/joined_kegg.txt
    
    export PROB_LIMIT=0.0000001
    export MAX_SIM_GENES=50
    for GENE_ID in `cut -f1 $E_DIR/vectors.txt`; do
        ##
        ## Дописываем к генам функции и пути
        ## gene_id<TAB>probability<function>
        ##
        
        join -t$'\t' \
            <(join -t$'\t' \
                <(cat ${E_DIR}/joined_kegg.txt \
                    | awk -F'\t' '{if (($2 == "'"$GENE_ID"'") && ($3<'"$PROB_LIMIT"')) print $1"\t"$3"\t"$4}' \
                    | sort -t$'\t' -k2,2 -g | head -n $MAX_SIM_GENES | sort -t$'\t' -k1,1 \
                ) \
                ${BASE_DATA_PATH}/5_pathways/gene_pathway_map.txt ) \
            <(cat ${BASE_DATA_PATH}/1_databases/genes_info.txt | cut -f1,4) \
        > ${E_DIR}/results_${GENE_ID}_limited_${PROB_LIMIT}_.txt
        
        ##
        ## probability<TAB>function
        ##
        #cat ${E_DIR}/results_${GENE_ID}.txt \
        #    | cut -f2,4 \
        #    | sort -t$'\t' -g -k1,1 -S100M \
        #    | sort -t$'\t' -k2,2 -u -s -S100M \
        #    | sort -t$'\t' -g -k1,1 -S100M \
        #    > ${E_DIR}/results_${GENE_ID}_functions.txt
        #cat ${E_DIR}/results_${GENE_ID}_functions.txt | head -n1000  > ${E_DIR}/results_${GENE_ID}_functions_top1000.txt
        #
        ##
        ## группируем гены по путям, сортируем группы по числу генов
        ##
        join -t$'\t' -2 3 \
            <(join -t$'\t' -1 2 \
                <(cat ${E_DIR}/results_${GENE_ID}_limited_${PROB_LIMIT}_.txt \
                    | cut -f4 | sort | uniq -c | awk '{print $1"\t"$2}' ) \
                <(cat ${BASE_DATA_PATH}/5_pathways/pathways_full.txt | cut -f1,2,4)
            ) \
            <(cat ${E_DIR}/results_${GENE_ID}_limited_${PROB_LIMIT}_.txt | cut -f2,3,4,5 | sort -t$'\t' -k3,3 -S100M) \
            | sort -t$'\t' -k2,2nr -k1,1 \
            | 6_annotation/group_genes_from_pathway.py \
            > ${E_DIR}/results_${GENE_ID}_pathways.txt
    done
    
    
    6_annotation/intersect_genes_avg \
        -m ${BASE_DATA_PATH}/5_pathways/pathway_vectors_avg_${E_VALUE}_.txt \
        -v $E_DIR/vectors.txt \
        -g $GENOMES_NUM \
        > $E_DIR/joined_pathways.txt

    join -t$'\t'  \
        <(cat $E_DIR/joined_pathways.txt | grep NC_000907_15301888225605723776 | cut -f2- | sort -t$'\t' -k1,1) \
        <(cat ~/data/5_pathways/pathways_full.txt | cut -f1,4) \
        | sort -k2,2 -g \
        | awk '{if ($2<0.0000001) print}' \
        | cut -f4 \
        | sort \
        | uniq -c | sort -nr -k1,1 \
        > $E_DIR/joined_pathways_uniq.txt
done
